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1 INTRODUCTION 

Granular materials are ubiquitous in nature and very common in industrial pro- 
cesses, but it is only recently that their unusual properties have begun to receive 
detailed attention from the physicists community ^ . The earliest documented 
studies of granular matter date back to Faraday ||, who discovered the convec- 
tive behavior of vibrated sand, and Reynolds ||, wno noted that compactified 
granular matter cannot undergo shear without increasing its volume. 

The behavior of vibrated granular matter in some aspects resembles that of 
a fluid, although there are crucial differences. Size segregation |J, for example, 
at first sight defies intuition. When a mixture of particles of different sizes is 
subject to vibration, the larger ones migrate to the top, irrespective of density. 
Also interesting is the layering instability |J of a binary mixture under pouring. 
Instead of a homogeneously mixed pile, under certain conditions an alternation 
of layers of both kinds of particles can be obtained. 

Similar demixing phenomena occur in granular materials subject to various 
kinds of external excitation. These seem to contradict the naive expectation that 
shaking should favor mixing, or take the system to a low-energy state. Many 
of the unusual properties of vibrated granular matter are in fact due to the 
dissipative character of interparticle collisions. An interesting example of the 
consequences of dissipation is inelastic clustering JjJ, by which particles tend to 
cluster together as their relative kinetic energy is completely lost during collisions. 
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The compactification of vibrated sand has been recently found to be logarith- 
mically slow ||, resembling glassy behavior, and a spin model with frustration 
has been proposed to model this process ||. This provides a bridge between the 
dynamics of spin glasses and vibrated granular matter. 

It is thus clear that granular materials present extremely interesting dynamic 
phenomena, but it is already at the much simpler level of static properties that 
unusual behaviors show up. Stress propagation in piles or packings of granulate 
matter has many uncommon features. When grains are held in a tall vertical 
silo, for example, the pressure at the bottom does not indefinitely increase with 
height but saturates after a certain value |10|. The excess weight is deviated 



towards the walls and equilibrated by friction forces. A related phenomenon is 
the formation of a pressure "dip" right below the apex of a conical pile of gran- 
ular matter 
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instead of the expected pressure maximum. These phenomena 
indicate that gravity-induced stresses do not propagate vertically but often de- 
viate laterally. Pressure saturation in silos and the pressure dip under piles are 
both due to the formation of "arches" O]. Many proposals to explain 

arching |T2|, [13] rest on the idea that friction plays an essential role, but recent 
studies [14, |l^, 16, 17] show that friction is not necessary. 

m ~ 



Photoelastic visualization experiments [18, | 



IT] show that stresses in 
granular matter concentrate along well defined paths. It is not clear whether 
the characteristic size of these patterns is finite, or limited by system size only. 
Stress-concentration paths are observable even on regular packings of monodis- 
perse particles |T^, [TjJ, their exact location being sensitively dependent on very 
weak particle irregularities. Stress-paths often suffer sudden rearrangement on a 
global scale when the load conditions are slightly changed [TBI EU1 E2| . For similar 



reasons, the fraction of the total weight that reaches the base of a silo can vary 
by large amounts under very weak perturbations, or when repeating the filling 
procedure with exactly the same amount of grain [13| . These phenomena demon- 
strate that slight perturbations can produce macroscopic internal rearrangements 
in granular matter. In other words, granular matter are internally unstable. 

In part because of the technological importance of the problem, and also be- 
cause of its interest from the point of view of basic science, much work has been 
done in recent years to understand the propagation of stresses in granular sys- 
tems. On the numerical side, several methods have been implemented. Classical 
Molecular Dynamics simulations [TBI, |2T], |23|. |24] , which usually include a fictitious 



damping term in order to allow the system to come to rest, are normally very 
cpu-intensive and thus limited to relatively small sizes. Alternatively, the elastic 
equations can be solved using symbolic software in order to obtain stress values 
which are free of numerical error [p^] . Lattice automata based on random con- 
tact disorder [26], are able to reproduce the observed dip under granular piles. 
Contact Dynamics simulations |22| provide an efficient way to include friction 



forces, and allowed numerical visualizations of stress concentration on relatively 
large systems. Lubricated Granular Dynamics jnj is a method to obtain the 
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equilibrium contact network of infinitely stiff networks and is based on the use of 
a fictitious damping with a singularity at zero distance. 

A large number of theoretical approaches to this problem are formulated on a 
, 129. and thus rest on the assumption that a length scale 



continuum 12, 27 



exists, below which fluctuations are negligible when compared to averages. It is 
not clear whether this assumption is easily justified for granular matter, where 
stress fluctuations seem to be at least as large as average stresses f2(| |2j] . 



A different type of modeling strategy starts by formulating a stochastic rule for 



stress propagation on a lattice, which is thereafter solved by various methods |21 
], p3|, or taken as the starting point for a continuum description 



24, BT 



In the simplest version of this approach [^T[ 24[ , only the vertical component of 
the transmitted force is considered, i.e. the problem is reduced to a scalar one. 
Despite the roughness of the approximation, this procedure gives good results for 
the average distribution of stresses P(w), in particular the observed exponential 
decay for large stresses 0, 21], 22[ . The occurrence of small stress is though 
strongly underestimated within this simple scalar model [21, |ST| . This is due to 
the fact that scalar "stresses" propagate vertically, with at most a diffusive width 
due to disorder, and therefore arching is not possible [31]. In order to correctly 
reproduce the small-stress part of P(w), which is arch-dominated, the vectorial 
nature of stresses has to be taken into account 



33|| . This brings in the 
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problem of stress signs, since now negative (traction) stresses, which do not exist 



on non-cohesive granular matter, cannot be easily avoided [31 



Figure 1: The contact network (right) associated to a granular pile is a graph in 
which nodes represent particles and are connected by an edge or bond whenever 
there is a nonzero (compressive) force between the corresponding particles. 

It is thus clear that granular matter has, from the point of view of its static 
properties, two noteworthy characteristics: 

• Stresses are not homogeneously distributed over the system but concentrate 
on paths that form a sparse network. 
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• The exact location of stress paths is susceptible of change under very weak 
perturbations, showing that granular matter is extremely unstable. 

Although there have been many proposals to describe stresses in granular 
matter, most of these models are largely phenomenological in nature and some- 
times contain unclear ad hoc assumptions. A deeper understanding of the above 
described particularities of granular matter has remained elusive. It is the pur- 
pose of this work to show that structural rigidity concepts can help us advance 
in this direction. 

We first review some structural rigidity notions in Section |||, and demonstrate 



in Section [2.1| that the contact network of a granular system becomes exactly 
isostatic in the limit of large stiffness. The consequences of this are discussed in 
Section |3|. An immediate consequence of isostaticity is the possibility of stress 
concentration, as briefly discussed in the beginning of Section ||. Most of the 
previous theoretical and numerical effort has concentrated on the description of 
stresses. There is though a complementary aspect to this problem, which has not 
been explored. This is the study of how displacements induced by a perturbation 
propagate in a granular system. We thus leave further discussion of stresses 
for future publications JT7|, and concentrate in understanding the behavior of 



induced displacements upon perturbation. This will lead us to the central results 
of this work, respectively: 

a) in Section [TT] it is shown that an isostatic phase transition takes place in 
the limit of infinite stiffness, and that the isostatic phase is characterized by an 
anomalously large susceptibility to perturbation. 

b) Section contains a discussion of the load-stress response function of the 
system in the light of a), which shows that isostaticity is responsible for the 
observed instability of granular matter. 

We will furthermore find that very large displacements are produced on iso- 
static networks when a site is perturbed. Section |3.3| clarifies the origin of these 
anomalously large displacements, while Section f| contains our conclusions. 



2 STRUCTURAL RIGIDITY AND GRANU- 
LAR NETWORKS 

The contact network of a frictionless packing of spherical particles can be defined 
in the following way (Fig. [I]): we let each particle center be represented by a 
point in space, and connect two of these points by a line (bond, or edge) when- 
ever there is a nonzero compression force between the corresponding particles. 
The networks so generated can be seen as particular cases of what is usually 
called frameworks in rigidity theory, i.e. structures made of points connected by 
rotatable bars. 

Structural rigidity studies the conditions that a network of points con- 
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nected by central forces has to fulfill in order to support applied loads, i.e. be 
rigid. The first studies of rigidity of structures from a topological point of view 



date back to Maxwell [3G]. Structural rigidity concepts were first introduced in 
the study of granular media by Guyon et al f37j , who stressed that granular con- 
tact networks differ from linear elastic networks in an important aspect: the first 
are only able to sustain compressive forces between grains. Technically speaking, 
force networks with a sign-constraint on stresses are called struts. Another typi- 
cal example of sign-constrained networks are spider webs, or cable structures, the 
elements of which (strings) can only sustain traction forces |38 |. Structures with 



interesting properties can be obtained by combining elements of both types, in 
which case the resulting framework is called a tensegrity structure [ 39f] . Several 
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Figure 2: a) A minimally rigid network in two dimensional space, composed 
of n points and 2n — 3 bars. If any bar in this network is removed, some of the 
points would cease to be rigidly connected to the rest. All bars in this network are 
therefore essential for rigidity, b) A redundant bar (dashed line) between points 
A and B can induce stresses ("self-stresses") in some part of the network (dark 
lines). The locus of self-stresses is the overconstrained subgraph associated with 
bar AB. This is exactly the set of bars which already provides a rigid connection 
between points A and B. If the length of bar AB is exactly equal to the distance 
d(A,B), i.e. if there is no length-mismatch, self-stresses will be zero. 

applications of rigidity related ideas and tools have been already presented in 
this book. Let us here only briefly refresh some concepts which we need for our 
discussion. 

A point in d dimensions has d degrees of freedom, while a rigid cluster has 
d(d + l)/2. Therefore if a set of n points forms a rigid cluster, it must be con- 
nected by at least dn — d(d + l)/2 bars. If a rigid cluster of n points has exactly 
dn — d(d + l)/2 bars, it is said to be (generically) isostatic, or minimally rigid 
(Fig. |2]a). A framework with more bars than necessary to be rigid is hyperstatic 
or overconstrained (Fig. ^b). Excess bars, which can be removed without intro- 
ducing new degrees of freedom, are called redundant. A bar is redundant when 
the two points it connects are already rigidly connected. Unless the length of the 
redundant bar is exactly equal to the distance between the points it connects, 
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self-stresses will be generated in some parts of the framework. Thus self-stresses 
are non-zero only within overconstrained subgraphs, and can be absent if there 
are no length-mismatches. 

We will discuss granular piles under the action of gravity, in which forces are 
transmitted to a supporting substrate. We could as well consider any other load 
condition, provided that an infinitesimal gravity field is added in order to remove 
indeterminacies in the positions of the particles. Because of the action of gravity, 
contact networks associated to static granular piles must be rigid. Otherwise 
some of the grains would be set in motion by gravitational forces. Because of the 
fact that grains are rigidly connected to a lower boundary, a redundant contact 
anywhere on the system will usually produce an overconstrained subgraph that 
extends all the way down to the rigid boundary. 
Some early attempts to numerically study the static properties of granular ma- 





Figure 3: In the limit of large stiffness-to-load ratio, i.e. when the compressive 
forces are small, or the rigidity large (left), the contact network of a granular 
packing is sparse and, as discussed in the text, becomes minimally rigid. If the 
compressive forces are increased, or the stiffness decreased (right), excess contacts 
(dashed lines) are created. 

terials have ignored the sign constraint, thus modeling granular piles as randomly 
diluted linear frameworks ffil . Due to this, it has been sometimes suggested that 



rigidity percolation concepts might be applied to granular networks f|T| . But 
this would require forces with power-law distribution, since the elastic percolation 



backbone is a fractal object at the critical point |4C], |42| , while experimental and 



theoretical studies |2l], |20], [22] suggest that force distributions display exponential 
decay on granular systems. 

This suggests that the sign-constraint associated with non-cohesive granular 
matter cannot be neglected. As we shall soon show, it is possible to see from a 
topological point of view that the sign-constraint has far-reaching consequences 
for the static behavior of granular aggregates. We demonstrate next that this 
restriction forces a stiff granular system to choose, among all possible equilibrated 
contact networks, only those with the specific topological property of minimal 
rigidity. 
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2.1 Isostaticity of stiff networks with a sign constraint 

Consider now a (i-dimensional frictionless granular pile in equilibrium under the 
action of external forces Fi (gravitational, etc) on its particles. We represent 
the contact network of this pile by means of a linear-elastic central-force network 
in which two sites are connected by a bond if and only if there is a nonzero 
compression force between the two corresponding particles in the pile (Fig. [[]). 

If the external compression acting on the pile is increased, particles will suf- 
fer a larger deformation, and therefore the number of interparticle contacts will 
increase (see Fig. |3|). Equivalently, if compression forces are weakened, or the 
stiffness of the particles made larger, the number of interparticle contacts, and 
thus the number of bonds on the equivalent contact network, will be reduced 
because there are no cohesive forces between particles. But there is a lower limit 
for the number of remaining contacts, given by the condition that all particles 
be rigidly connected, otherwise they would move until new contacts are estab- 
lished. Therefore one may expect that in the limit of infinite stiffness the resulting 
contact network will be minimally rigid. Let us now formalize this observation. 

Because of linearity, stresses on the bonds of the linear-elastic equivalent 
network can be uniquely decomposed as 

h iff ■ f!r d a) 

, where /^ elf are self-stresses, and / i 1 ° ad are load-dependent stresses. 

Load-dependent stresses are linear in the applied load: if all loads are rescaled, 
fl° ad are rescaled by the same factor. But fl° ad are not changed if all elastic 
constants, or stiffnesses, are multiplied by a factor. 

Self stresses in turn do not depend on the applied load. They are in general 
linear combinations of terms of the form kijSij, where kij is the stiffness of bond 
ij and its length- mismatch (See Fig. ^j. The length mismatch of a bond is 
defined as the difference between its repose length and its length in an equilibrium 
configuration under zero external loads. 

As already discussed, length mismatches, and thus self-stresses, only arise within 
overconstrained subgraphs |34|, |3^]: those with more contacts, or bonds, than 



strictly necessary to be rigid. For if a graph is not overconstrained, then all its 
bonds, regarded as linear springs, can simultaneously attain their repose lengths 
while still being in equilibrium under zero external load. 

It is easy to see that a bounded overconstrained subgraph with nonzero self- 
stresses must have at least one bond subject to a negative (traction) stress: As 
discussed, self-stresses are equilibrated without the intervention of external forces. 
It then suffices to consider a joint belonging to the envelope of the overconstrained 
subgraph. Since bonds can only reach it from one side of the frontier, self-stresses 
of both signs must necessarily exist in order for this joint to be equilibrated. 

Now imagine rescaling all stiffnesses according to k — > Xk (both in the granular 
pile and in our equivalent elastic system). In doing so, all self-stresses are rescaled 
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Figure 4: a) In order for seven particles to be in contact forming an overcon- 
strained graph as shown here, one condition has to be satisfied by the radii. One 
is free to choose the radii of 6 of them, but the seventh will be uniquely deter- 
mined by this choice, b If one of the particles is for example slightly too large to 
satisfy the required condition, one contact will be opened, restoring isostaticity. 
In this example the contact between particles 2 and 7 is open, but any other bond 
between the central particle and its neighbors could have been chosen. If on the 
other hand particle 2 were too small, one of the contacts between the external 
particles would be lost. 

by A, but the load-dependent stresses remain constant. Therefore, if self-stresses 
were non-zero, in the limit A — > oo at least one bond of the network would have 
negative total stress (eq. ([I])). This is not possible since traction forces are not 
allowed on our granular pile by hypothesis. Thus the existence of self-stresses 
is not possible in the limit of large stiffness if there is a sign-constraint on total 
stresses. 

For reasons already reviewed, in order for self-stresses to be zero, one or more 
of the following conditions will have to be satisfied by granular contact networks 
in the limit of large stiffness: 

a) to have all length mismatches equal zero on overconstrained graphs. 

b) to have no overconstrained graphs at all (the network is isostatic). 

Condition a) requires that, even when overconstrained graphs exist, particle 
radii satisfy certain conditions in order to exactly fit in the holes left by their 
neighbors (see Fig. [|). But as soon as particles have imperfections or polydis- 
persity (no matter how small if their rigidity is large enough) self-stresses would 
appear if the contact network is hyperstatic. In other words condition a) cannot 
be generically satisfied, if by generic we understand for a "randomly chosen" 
set of radii. Therefore condition b) must generically hold, i.e. there will be no 
overconstrained subgraphs in the limit of large stiffness-to-load ratio. 

From the point of view of this work, most experimentally realizable pack- 
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ings fall under the category "generic", since small imperfections in radius are 
unavoidable in practice. 

In view of the above, we can now conclude that: 

the contact network of a polydisperse granular pile becomes isostatic 
when the stiffness is so large that the typical self-stress, which is of 
order ke would be much larger than the typical load-induced stress. 

Exceptions to this rule are packings with periodic boundary conditions, be- 
cause they are not bounded, and packings in which the radii satisfy exact condi- 
tions in order to have zero length-mismatches in overconstrained graphs, because 
they have no self-stresses E| . One can for example consider a regular triangular 



packing of exactly monodisperse particles, in which case the associated contact 
network will be the full triangular lattice |25[], i.e. hyperstatic. But the contact 



network (and the properties of the system as we shall soon see) will be drastically 
modified as soon as a slight polydispersity is present [F23|, [TBI, [H3[] if the stiffness 
is large enough. Therefore, while one is free to consider specific packings which 
are not isostatic, in practice these cannot be realized for hard particles, since 
any amount of polydispersity, no matter how small, will force some contacts to 
be opened so as to have an isostatic contact network. We now discuss how this 
affects the properties of granular systems. 



3 CONSEQUENCES OF ISOSTATICITY 

Isostaticity has been sometimes imposed in numerical models [26|, as a condi- 
tion allowing one to calculate stresses by simple propagation of forces. Recently 



isostaticity was reported by Ouaguenouni and Roux [14]], who use an iterative 
numerical algorithm to find the stable contact network of a set of rigid disks. Our 
discussion in the previous section shows that isostaticity is a generic property of 
stiff packings, and appears because negative stresses are forbidden (an equivalent 
conclusion would be reached if only traction stresses were allowed f38|). We will 
now show that isostaticity has important consequences on the way the system re- 
acts when it is perturbed, but before starting a more rigorous analysis, let us first 
discuss some of the most important differences between isostatic and hyperstatic 
systems, on an intuitive level. 

Imagine perturbing an elastic network by letting an equilibrated pair of collinear 
forces act on the ends of a given bond, and consider how stresses and equilib- 
rium positions are modified in the whole system. A properly chosen change in 
the repose length of this bond would have exactly the same effect, so we can 
alternatively consider the perturbation to be a change in length or a couple of 
forces. 

• On overconstrained rigid networks, stresses are correlated over long dis- 
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tances: if a bond is stretched as described above, stresses will be modified 
on other bonds far away from the perturbation. This is so because self- 
stresses percolate through the system. 

But for the same reason the displacements of the sites from their original 
positions, induced by this local perturbation, decay very fast with distance. 
The reason for this is that self-stresses oppose the perturbation and thus 
tend to "quench" its effect. 

• On isostatic systems, stresses are uncorrelated. If we change the length of 
(or apply a equilibrated pair of forces to) an arbitrary bond, stresses on all 
other bonds remain unchanged, because they only depend on external loads, 
and not on bond lengths. This is a trivial property of isostatic systems. 

But the equilibrium position of many sites will be in general modified if 
one of the lengths is changed, and therefore displacements induced by a 
perturbation may be felt far away from its origin. 



Figure 5: Part of the stress-carrying backbone for rigidity percolation on a ran- 
domly diluted triangular lattice. Overconstrained rigid clusters (thick lines) are 
isostatically connected to each other by cutting bonds (thin lines). Cutting bonds 
are critical in the sense that they provide a minimally rigid connection. All of 
them are essential for rigidity. 

Therefore on a hyperstatic system, perturbations of stresses propagate over 
long distances, while on isostatic systems it is the displacement field which may 
display long-range correlations. Stresses are uncorrelated on isostatic systems, 
and thus arbitrarily large stress gradients are possible. Hyperstatic systems, on 
the opposite hand, have smoothly varying stresses because of the strong corre- 
lations introduced by overconstraints. This provides some indication that stress 
concentration is possible because of isostaticity. 
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We can get further insight into the meaning and potential consequences of iso- 
staticity from recent studies of central-force rigidity percolation |42|. Rigid back- 



bones, the stress-carrying components of rigidly connected clusters, are found 
to be composed of large overconstrained clusters, is o statically connected to each 
other by critical bonds (also called red bonds, or cutting bonds - See Fig. H). 
Overconstrained clusters have more bonds than necessary to be rigid, so any one 
of them can be removed without compromising the stability of the system. But 
the rigid connection among these clusters, provided by critical bonds, is isostatic, 
or minimally rigid. In other words, cutting one critical bond is enough to produce 
the collapse of the entire system, because each isostatic bond is by definition es- 
sential for rigidity. Thus we may expect that stretching a critical bond will have 
a measurable effect on a large number of sites. 

But in percolation backbones, critical bonds only exist very close to the rigid- 
ity percolation density p c . Above p c there is percolation of self-stresses [[£2| and 
thus the rigid backbone is hyperstatic. Even exactly at p c the number of critical 
bonds is not extensive, but scales as L 1 ' v where v is the correlation-length expo- 
nent [42, [HJ]. Consequently critical bonds are relatively few at p c , and virtually 



absent far from p c . Thus if we perturb (cut or stretch) a randomly chosen bond 
in a percolation backbone, most of the times the effect will be only be local since 
no critical bond will be hit. 

The important new element in stiff granular contact-networks is the fact that 
all contacts are isostatic, i.e. there is extended isostaticity. In this case, if any of 
the bonds (contacts) where removed, the pile would cease to be rigidly connected 
to the supporting boundary below it. Because of this, if the length of any of 
the network's bonds is changed (which corresponds to a variation in one of the 
particles radii) the equilibrium position of a finite fraction of the particles will 
also be changed. For these reasons, one may expect that isostaticity will produce 
a large sensitivity to perturbation in granular networks. 



3.1 Susceptibility to perturbation 

We now quantify the degree of susceptibility to perturbation, and then see whether 
the intuitively appealing ideas we have just discussed are in fact verified on spe- 
cific models. In order to provide a formal definition of susceptibility, we introduce 
an infinitesimal change in the length 1^ of a randomly chosen bond of the net- 
work, and record the induced displacement 5i suffered by all particle centers. We 
then define the system's susceptibility D as 

N 

£ = E? (2) 

1=1 

, where N is the total number of particles on the system. These measurements are 
done for variable amounts of over constraints (excess contacts) randomly located 
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Figure 6: A linear chain of springs in which overconstraints (diagonals) are 
present with probability O v is useful as a toy model to understand the influence 
of isostaticity on the propagation of perturbation. The characteristic distance 
for the decay of displacements induced by a perturbation (a bond-stretching) 
diverges in the limit O v — > and thus the susceptibility to perturbation also 
diverges (see text). 



on an elastic network, and averages are performed over disorder. In this way 
D(O v ) is obtained, where O v is the density of overconstraints. According to our 
previous discussions, we expect D to increase as O v — > 0, which is the isostatic 
limit. 

We start by discussing a toy model shown in Fig. || a quasi one-dimensional 
system composed of linear elastic bonds, in which diagonals are present with 
probability O v . The system with no diagonals (O v = 0) is exactly isostatic, 
therefore each diagonal is an overconstraint, or redundant bond. We assume for 
simplicity that all bonds have the same stiffness k and length I (diagonals have a 
length y/2l). We are interested in calculating the average horizontal displacements 
u{x) induced by a length change A/ in the left-most horizontal bond, as a function 



of the density O v of overconstraints (diagonals). A simple calculation |17| shows 
that, after averaging over disorder, the displacement field u(x) satisfies 

d 2 u 

— = #O v u (3) 

, where k is some constant. Therefore u(x) = u(0) exp{— k0^ 2 x} and we see that 
there is exponential decay with distance, with a characteristic length £(O v ) ~ 
0~ l l 2 . This "persistence length" diverges at O v = 0, which corresponds to the 
isostatic limit. Consequently O v = is a critical point, and D as defined above 
is divergent there. The divergence of D in this model is linear with system size. 
Now let us see whether isostaticity has comparable effects in two dimensions. In 
the spirit of previously proposed models |26j we consider a triangular packing 
oriented as in Fig.^a, made of very stiff disks with small polydispersity, under 
the action of gravitational forces. The polydispersity is assumed to be small 
enough such that disk centers are approximately located on the sites of a regular 
triangular lattice, and the stiffness to load ratio large enough such that the contact 
network is isostatic, according to our discussion in the previous section. 
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Figure 7: a) In two dimensions, a triangular packing is used in order to nu- 
merically measure susceptibility to perturbation. The two first layers of particles 
(shaded) is regarded as a fixed rigid boundary supporting the load of the upper 
ones. If a site is shifted, only particles within its "cone of influence" (dashed line) 
can be displaced. In this example, 6 layers of 6 particles each are displayed, b) 
Appropriately choosing among these three isostatic configurations for each site, 
only compressive stresses are produced. First S is chosen with probability 1/2. 
If S is not chosen then either R or L are, depending on the sign of the horizontal 
force acting on the site (see text). 

The full problem of generating realistic contact networks that respect the 
constraint of no traction forces is a difficult one. Several approaches have been 
proposed f41| , Pq| , all involving some degree of approximation even for geomet- 
rically simple settings. In all these models, the triangular lattice was oriented 
with one of its principal axis horizontal, i.e. normal to gravity. There is though 
some advantage in considering a different orientation, such that one of the prin- 
cipal axis of the lattice is parallel to gravity (Fig. ^). In this case there is no 
need for recursive checks of positiveness of stresses since the disordered contact 
network can be built in a fashion that guarantees positive stresses. Our model 
is defined as follows: We ask that each site be supported by exactly two out of 
its three lower neighbors, thereby ensuring that only isostatic contact networks 
are generated. This condition gives three possible local configurations which we 
call left (L) symmetric (S) and right (R) respectively and are depicted in Fig.[7|b. 
Choosing a local configuration of bonds on each site produces a sample, or one 
realization of the disorder. Clearly not any choice of bond configurations give 
rise to a contact network with positive stresses. But it is possible to satisfy the 
positivity constraint and still have disorder in the following way: 

1. For each site, starting from the uppermost layers and proceeding down- 
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wards, we choose configuration S with probability 1/2 

2. If S was not chosen, then either R or L are, according to the sign of the 
horizontal force- component F x acting on that particle: if F x > (F points 
rightwards), R is chosen. If F x < 0, then L is chosen. If the horizontal 
component is zero then R or L are chosen at random. 

Our model has no geometrical disorder, and this is justified by our assumption 
of small polydispersity, but we keep contact disorder and isostaticity which are the 
important characteristics of real granular networks in the limit of large stiffness. 

There is no reason to think that the method we have chosen generates all 
possible equilibrated contact networks that satisfy isostaticity and positiveness 
of stresses. It seems in principle possible to have some sites making contact to all 
three downward neighbors and still have isostaticity, by simultaneously opening 
some other contacts. But our aim here is not to provide a realistic model for 
granular contact networks but to test whether isostaticity has important effects 
on the properties of a two-dimensional network. 

In order to accomplish this we will measure, on the networks so generated, the 
susceptibility defined above, and compare the results with those obtained on sys- 
tems with a finite density O v of overconstraints randomly located on the network. 
A non-zero density of overconstraints O v mimics, as discussed previously, the ef- 
fect of increasing the mean pressure on the system (or reducing the stiffness), 
since this would produce a larger number of contacts, in excess of isostaticity, to 
be established between particles. 

A finite density O v of overconstraints is introduced in this model by letting all 
three bonds be connected below a given site, with probability O v . Each third bond 
introduced in this way creates an overconstrained subgraph that extends all the 
way down to the rigid boundary. The limit O v = 1 gives the fully connected trian- 
gular lattice, which of course has no disorder. After building a contact network 
with a specified density of overconstraints as described above, an infinitesimal 
upwards shift is introduced in a randomly chosen site on the lowest layer, and 
the induced displacement field Si is measured. 



If the network is isostatic (O v = 0) one can calculate all stresses |4T], ^6 



and displacements |T7[ in a numerically exact fashion so that systems of 2000 x 
2000 particles may be simulated on a workstation. The idea is that stresses 
are propagated downwards and displacements upwards. The way in which the 
induced displacements are propagated upwards is easily calculated [|I7| by noting 
that, when the network is isostatic, all bond lengths (except the perturbed one) 
must remain constant. On the other hand when the network is overconstrained 
(O v > ), stresses and displacements can no longer be exactly calculated. In this 
case one has to solve the elastic equations in order to find the new equilibrium 
positions after the perturbation. This is done in the limit of linear elasticity 
(since the perturbation is infinitesimally small) by means of a conjugate gradient 
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Figure 8: a) Susceptibility D y (defined in the text), versus density of overcon- 
straints O v , as measured on two-dimensional triangular packings of total height 
H = 12, 24, 50 and 100 layers. O v = is the isostatic limit, and corresponds to 
granular packings with infinite stiffness, as demonstrated in the text, b) Sus- 
ceptibility D y versus system height H, for fixed fractions O v of overconstraints: 
0.00 (empty circles), 0.01 (squares), 0.02 (diamonds), 0.05 (triangles), 1.00 (full 
circles). It is clear in this figure that, for O v > , D y saturates to a finite value 
in the H — > oo limit. For isostatic systems, on the other hand, D y diverges 
exponentially fast with H . See also Fig ^|a). 



solver. In this case the calculations are much more time-consuming so that only 
relatively small systems, or order 100 x 100 particles can be studied if O v ^ 0. 
Supercomputers are required for this part of the calculation [ 4"S| . A cross-check 
of the computer programs was done by comparing the results obtained with the 
direct solver for isostatic systems, with those produced by the conjugate gradient 
solver with no overconstraints, on systems of up to 100 x 100 particles. Excellent 
agreement was found in all cases, for stresses as well as as for displacements. 

In this way the susceptibility D y (H,O v ) =< Y^i^Ui 2 > is measured, where 
Syi is the vertical displacement of site i due to the perturbation, and <> stands 
for average over disorder realizations. The system consists of H layers of H grains 
each, so that N = H x H . Figure ||a shows the susceptibility D y as a function of 
the density of overconstraints O v , for several system heights H. We see that D y 
increases rapidly on approach to the isostatic limit O v = 0, and that this increase 
is faster for larger systems, meaning that D y diverges at O v = in the H — » oo 
limit. Figure ||b shows the same data now plotted as a function of system size, 
for several values of the density of overconstraints. For any O v ^ 0, D y goes to a 
finite limit for large sizes, while it diverges with system size if O v = 0. Data for 
much larger systems can be obtained in the isostatic case using the direct solver 
program, and are displayed in Fig. |9|a. This plot shows that D y {O v = 0) is of the 
form log D y ~ H, that is, D y diverges exponentially fast with system size. 

These numerical results demonstrate that a phase transition occurs at O v = 0, 
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Figure 9: a) D y grows exponentially with height when O v = 0. b) The prob- 
ability Ph(8y) to have an induced displacement 8y is power-law distributed on 
isostatic networks. Results are shown at h — 200, 400, 600, . . . , 2000 layers above 
the perturbation, for positive values of Sy only. Ph(5y) is approximately symmet- 
ric, and has a finite peak of weight 0.47 at 5y = 0. Thus at any height above the 
perturbation, approximately one half of the sites within the influence cone are 
not vertically shifted. 

where anomalously large susceptibility sets in. In a way which is consistent with 
our intuitive expectations in the previous section, and with the one dimensional 
toy model, isostaticity is also in two dimensions responsible for a large suscepti- 
bility to perturbation. An important and surprising difference is the fact that, 
at the isostatic critical point O v = 0, the susceptibility D increases exponentially 
fast with system size, whereas it grows only linearly with size in one dimension. 

Seeking to understand the surprisingly fast growth of D with system size, 
the probability distribution Ph(o~y) to have a vertical displacement 5y, h layers 
above the perturbation has been measured on isostatic systems of 2000 x 2000 
particles. Only sites within a 120 degree cone whose apex corresponds to the 
perturbed bond may feel the effect of the perturbation (Fig. ^a). Ph(5y) thus 
gives the probability for a randomly chosen site inside the influence cone of the 
perturbed bond and h layers above it, to have a vertical displacement 5y. Sites 
outside this cone have 5 = 0, so our measurements essentially correspond to a 
system of infinite width. 

Figure |b shows the result of these measurements at h = 200, 400, 600, . . . , 2000 
layers above the perturbation. Only positive values of S are displayed in this fig- 
ure since Ph(8y) is approximately symmetric. Ph(8y) is found to be consistent 
with a power-law behavior with an /i-dependent cutoff: 

P h {8y) ~ h-r\5y\- e (4) 

for 5y < 8 M {h). 

It is also evident from Fig. ^|d, that the cutoff 5m (h) grows exponentially with 
increasing distance h from the perturbation. Fitting the curve corresponding to 
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h = 2000 in the interval 10 10 < 5y < 10 30 , an estimate 6 = 0.98 is obtained, 
suggesting that — ► 1 asymptotically. Normalization then requires p = 1, since 
the cutoff 5m (h) increases exponentially with h. 

Similar measurements where done for (smaller) systems with a finite density 
of overconstraints O v , in which case the distribution of displacements presents a 
size-independent bound (Fig. |TT|b). 

Thus, in two dimensions an isostatic phase transition takes place at O v = 0, 
and the resulting isostatic phase is characterized by a susceptibility to perturba- 
tion that grows exponentially fast with system size. The distribution of induced 
displacements is power law, with a cutoff that grows exponentially with distance 
from the perturbation. Of course one does not expect to be able to really measure 
exponentially large values of displacements on granular systems. The calculations 
reported here are valid for infinitesimal displacements, and in with this in mind 
the contact network is considered to remain unchanged during the perturbation. 
In practice, internal rearrangements would occur before we could detect very 
large values of displacement on a real pile. So how can we know if the huge 
susceptibility to perturbation that our calculations predict have any observable 
effect? This is is discussed in the next section. 



3.2 Isostaticity implies instability 



In order to clarify the relevance of the findings described in section |3.1| in relation 
with the observed unstable character of granular packings [0, [19], ^0J, we 
must first demonstrate the equivalence between induced displacements on site % 
and the load-stress response function b) of the stretched bond b with respect 
to a load on site i. 

The network's total energy can be written as 

N 

E = Y J W l y % + 1/2 £ h(l b ~l° b ) 2 (5) 

i=l b 

, where the first term is the potential energy (Wj are particle's weights) and the 
second one is a sum over all bonds and accounts for the elastic energy. lb are the 
bond lengths in equilibrium and 1® their repose lengths (under zero force). Upon 
infinitesimally stretching bond b', equilibrium requires that 

E^# + £U^-O§^ = o (6) 

i OLy m Ol b i 

where the second sum goes over bonds ov that belong to the same overconstrained 
graph as b' does. This is so since bonds not overconstrained with respect to b' do 
not change their lengths as a result of stretching b'. Since stress /& on bond b is 
fb = kb(l b — lb) this may be rewritten as 
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If there are no overconstrained graphs the left hand sum only contains bond b 
itself, therefore, 

A = 1^1 (8) 

showing that, in the isostatic case, the displacement 8yf' = |p induced on site i 
by a stretching of bond b is the response function Q(i, b) of stress with respect 
to an overload on site i. 

Taking averages with respect to disorder on equation (||), we obtain 

<f b >=J2W l <6yl b) > (9) 

i 

, and since average stresses on a given layer grow linearly with depth, we must 
have 

< SyV > H ~ H- 1 (10) 

We have seen that the second moment of Pn(8y) diverges as exp{H}, while ([TOD 
shows that its first moment goes to zero with increasing H . This can only hap- 
pen if Pn(8y) is approximately symmetric (this is numerically verified), which 
demonstrates that large positive and negative values of 5y appear with similar 
probability. Given now the equivalence between induced displacements and the 
load-stress response function, the existence of large negative induced displace- 
ments means that a positive overload at a random site i, would often produce a 
(very large) negative stress on any arbitrarily chosen bond b. This in turn indi- 
cates that the system will have to rearrange itself in order to restore compressive 
forces, since negative stresses are not possible. In other words, isostatic pack- 
ings are unstable to small perturbations, and will reorganize themselves on the 
slightest change in load, in order to find a new stable (compression only) contact 
network. 

In order to finish the demonstration that instability is a consequence of iso- 
staticity, we still have to show that a finite density of overconstraints would make 
the response function bounded again. 

When there are overconstraints, Pn{8y) is no longer critical but bounded as 
our numerical simulations show (see Fig. [II]). But in this case by is no longer 
identically equal to the load-stress response function, i.e. ([8]) no longer holds. One 
can nevertheless see, by looking at formula (0), that the weight-stress response 
function Q(i,b) must be bounded if Pn(8y) is. Therefore in the overconstrained 
case a finite overload of order < / > is necessary in order to produce rearrange- 
ments, and the system is thus no longer unstable. 

3.3 Pantographs 

The exponential growth of 5m{H) is responsible for the observed exponential 
behavior of the total susceptibility in the isostatic case. But we have yet to 
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understand for which reason exponentially large values of displacements do ex- 
ist. Surprisingly this can be explained in very simple terms. The appearance 




Figure 10: The observed exponential growth of induced displacements with 
distance to the perturbation is due to the existence of "pantographs" as the one 
shown in this figure. Upon stretching bond B by a small amount 5, site A moves 
vertically by an amount 25. Conversely a unitary weight at A produces a stress 
of magnitude 2 on bond B. This multiplicative effect only exists on isostatic 
systems, and is lost if the network is overconstrained. In the overconstrained 
case, a stretching of bond B would generate internal stresses on other bonds that 
oppose the deformation, and the displacement of site A would be much smaller. 



of exponentially large values of displacements is due to the existence of "lever 
configurations" or "pantographs" , which amplify displacements. 

Fig. ITDI shows an example of a pantograph with amplification factor 2. When 
the dashed bond is stretched by e, site A is vertically shifted by 2e. Given that 
there is a finite density of similar pantographs on the system, it is clear that 
displacements will grow exponentially with system height. 

This amplification effect only exists in the isostatic limit: Pantographs as the 
one in Fig. |10| are no longer effective if blocked by overconstraints. For exam- 
ple, an additional (redundant) bond between site A and the site below it would 
"block" the amplification effect of the pantograph. Then, a unitary stretching 
of bond B would induce stresses in the whole pantograph, but only a small dis- 
placement of site A. 
In order to understand why the transition occurs at zero density of overcon- 
straints and not at any finite density, it is extremely important to notice that 
pantographs are composed of all sites suffering displacement when the perturbed 
bond is stretched. Thus a typical pantograph covers a finite fraction of the sys- 
tem, and any non-zero fraction of redundant bonds is enough to place at least one 
excess bond on it, eliminating the lever effect. This explains why anomalously 
large induced displacements only exist in the isostatic limit O v = 0. 
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Figure 11: The effect of isostaticity is dramatically illustrated by a comparison 
of Ph{Sy) with and without overconstraints. Here by the displacement of a site, 
induced by a bond- length perturbation h layers below it. For any nonzero density 
of overconstraints (left, O v = 0.02 in this case), induced displacements decrease 
with distance. This is the usual behavior on an elastic continuum. On the 
other hand if O v = 0, i.e. when the system is isostatic (right), the distribution 
of induced displacements gets broader when the distance h to the perturbation 
increases. This is due to the multiplicative effect of pantographs (see Fig. p~0|) . 
These results where obtained on systems of total height 100 layers. 

4 CONCLUSIONS 



We have shown that the contact network of granular packings becomes exactly 
isostatic in the limit of large stiffness-to-load ratio, i.e. when the stiffness is large 
or the mean compressive load is small. We have furthermore provided analyti- 
cal (in Id) and numerical (in 2d) evidence that isostaticity is responsible for the 
appearance of a large susceptibility to perturbation, defined as the sum of the 
square displacements induced by a small bond-stretching. When an arbitrary 
bond is stretched on an overconstrained system, the effect of this perturbation is 
only felt locally. On the contrary, on an isostatic system the induced displace- 
ments grow with distance. This surprising phenomenon makes the susceptibility 
diverge exponentially fast with system size, and is produced by the existence of 
"pantographs" : network mechanisms that amplify displacements in the same way 
a lever does. 

We have also clarified the relationship between the susceptibility to perturba- 
tion defined in this work and the experimentally observed instability of granular 
networks. This was done using an equivalence between induced displacements 
and the weight-stress response function. The existence of negative values for the 
response function and the relation of this fact with instability were first discussed 
in the context of a phenomenological model for stress propagation JJT]. In that 
model, the appearance of large negative values for the response function (cor- 
rectly identified as a signature of instability by the authors) is a consequence of 
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ad-hoc assumptions about the way in which stresses propagate downwards. That 
work though does not correctly identify the physical origin of instabilities. The 
vectorial character of the transmitted quantity is not the reason (while it is a 
necessary condition), as easily illustrated by an overconstrained network. The 
reason by which granular contact networks are unstable is that they are isostatic. 

Thus stiffness produces isostaticity. Isostaticity is responsible for the tendency 
to global rearrangement upon slight perturbation of stiff granular materials. Any 
non-zero density of overconstraints is enough to destroy criticality and therefore 
drastically reduce instabilities. Therefore "soft" granular packings are stable. 
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